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We present a Bayesian model selection analysis of WMAP3 data using our code CosmoNest. We 
focus on the density perturbation spectral index ns and the tensor-to-scalar ratio r, which define 
the plane of slow-roll inflationary models. We find that while the Bayesian evidence supports the 
conclusion that ns 7^ 1, the data are not yet powerful enough to do so at a strong or decisive level. 
If tensors are assumed absent, the current odds are approximately 8 to 1 in favour of ns 7^ 1 under 
our assumptions, when WMAP3 data is used together with external data sets. WMAP3 data on its 
own is unable to distinguish between the two models. Further, inclusion of r as a parameter weakens 
the conclusion against the Harrison-Zel'dovich case (ns = 1, r = 0), albeit in a prior-dependent 
way. In appendices we describe the CosmoNest code in detail, noting its ability to supply posterior 
samples as well as to accurately compute the Bayesian evidence. We make a first public release of 
CosmoNest, now available at www.cosmonest.org. 



PACS numbers: 98.80.-k 

I. INTRODUCTION 

The recent three- year results from WMAP pj have pro- 
vided the first firm indications that the spectral index 
of primordial density perturbations, ns, differs from the 
Harrison-Zel'dovich case ns = 1- The likelihood func- 
tion in models with varying ns suggests that ns = 1 is 
excluded at around three to four sigma, in cosmologies 
with no significant tensor contribution to the microwave 
anisotropies. However, the WMAP team stress that their 
result, based on a chi-squared per degrees of freedom ar- 
gument, needs to be checked using the more sophisticated 
technique of Bayesian model selection Q H- That 
is the aim of the present paper, building on our previ- 
ous analysis of the WMAP first-year data using our code 
CosmoNest |5j. 

We will consider two different scenarios. The first 
concerns the spectral index alone, under the assump- 
tion that there are no primordial gravitational waves 
(parametrized by the tensor-to-scalar ratio r) . As there is 
presently no indication for gravitational waves, this anal- 
ysis addresses the question of whether ns 7^ 1 should be 
considered part of the standard cosmological model. Sec- 
ondly, we consider the plane of slow-roll inflation models 
parametrized by ns and r, representing the simplest class 
of inflation models (for an extensive review, see Ref. |||). 
This latter analysis determines the extent to which slow- 
roll inflation models have benefitted from the new data. 



II. BAYESIAN MODEL SELECTION 

One of the most important classes of statistical prob- 
lems in science, and particularly in cosmology, is deter- 
mining the best fit to data in the case where the un- 
derlying model (i.e. the set of parameters to be varied) 
is unknown. Typically, each parameter represents some 
physical effect that may influence the data, but one can- 
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not simply include all possible physical effects simulta- 
neously, as the data may be insufficiently constraining 
and all parameters become undetermined and biases get 
introduced Q. In the Bayesian framework, the solution 
is model selection statistics, which set up a tension be- 
tween goodness of fit to the data and model complexity. 
A model selection statistic does not care about the pre- 
ferred values of the parameters defining the model, but is 
a property of the model itself, where here 'model' means 
both a choice of the set of parameters to be varied and 
the prior ranges for those parameters. 

A key application of model selection is to provide a 
robust criterion for judging when data requires the ad- 
dition of new parameters. Many of the most pressing 
questions in contemporary cosmology are of this type, 
such as whether the dark energy density evolves with 
redshift, or whether primordial gravitational waves exist. 
For the present data compilation following the WMAP3 
announcement, it is the spectral index which is placed in 
the most interesting position — does WMAP3 convinc- 
ingly exclude the possibility that ns is precisely unity, as 
conjectured by Harrison and by Zel'dovich || long before 
the inflationary mechanism was discovered? In model se- 
lection terms, does the improved fit that a varying ns 
allows justify its inclusion as an extra variable fit param- 
eter? 

One might wonder why we should bother with a model 
selection analysis of a result which a parameter estima- 
tion analysis says is already at three to four sigma level. 
The answer is that this significance level is exactly where 
model selection techniques are at their most crucial. It 
has long been recognized in the statistics community that 
Bayesian model selection analyses can give results in con- 
tradiction with inferences based on 'number of sigmas'; 
this is known as Lindley's 'paradox' and is nicely sum- 
marized by Trotta [10(. Basically, Bayesian inference is 
inconsistent with the idea that there is a universal thresh- 
old, such as 95%, beyond which results should be seen as 



2 



definitive; instead such a threshold should depend both 
on the data properties and the prior parameter ranges of 
the models being compared. The Lindley paradox usu- 
ally manifests itself for results with significance in the 
range two to four sigma |10| . which as it happens is ex- 
actly where WMAP3 has placed ns- 

In a full implementation of Baycsian inference, the key 
statistic is the Bayesian evidence E (also known as the 
marginalized likelihood), which has the literal interpre- 
tation of the probability of the data given the model 
0, 0, H ■ According to Bayes theorem, it therefore up- 
dates the prior model probability to the posterior model 
probability. It is simply the average of the likelihood over 
the prior parameter space. Often, the quantity of inter- 
est is the ratio of evidences of two models Mi and Mo, 
called the Bayes factor and denoted Bio, which indicates 
how the relative model probabilities have been updated 
by the data. The evidence has been exploited in a range 
of cosmological studies [jj CU Ell • 

Computing the evidence is more challenging than cal- 
culating parameter uncertainties, as it requires knowl- 
edge of the likelihood throughout the prior parameter 
volume rather than only in the vicinity of its peak. So 
far brute force methods such as thermodynamic integra- 
tion, though accurate, have proved to be computationally 
very intensive |l2j . while approximate information crite- 
rion based methods often lead to results which do not 
agree and hence can be ambiguous 0, El • We have re- 
cently developed an implementation of an algorithm due 
to Skilling known as Nested Sampling |l4j , which we call 
CosmoNest Q, which is able to carry out such calcula- 
tions efficiently. It is a Monte Carlo method, but not a 
Markov chain one. We describe the code extensively in 
Ref. [j| and in the appendices of this article. 

In assessing the significance of a model comparison, 
a useful guide is given by the Jeffreys' scale Q- La- 
belling as Mi the model with the higher evidence, it 
rates hi -Bio < 1 as 'not worth more than a bare men- 
tion', 1 < Initio < 2.5 as 'substantial', 2.5 < ln£>io < 5 
'strong' to 'very strong' and 5 < In B\q as 'decisive'. Note 
that InSio = 5 corresponds to odds of 1 in about 150, 
and ln£>io = 2.5 to odds of 1 in 13. 



III. APPLICATION TO WMAP3 

Throughout we use a data compilation of the WMAP3 
TT, TE and EE anisotropy power spectrum data 0, 
together with higher I CMB temperature power spec- 
trum data from ACBAR [3, CBI [H, VSA P3, and 
Boomerang 2003 and also matter power spectrum 
data from SDSS H and 2dFGRS H3- Following the 
approach of Ref. |2l|. we use the updated beam error 
module, and do not marginalize over the amplitude of 
SZ fluctuations. For the higher I CMB data, we neglect 
those bands that overlap in I range with WMAP (as in 
Ref. Q), so that they can be treated as independent mea- 
surements. 
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FIG. 1: Marginalized likelihood of n s for WMAP alone 
(dashed) and WMAP+all (solid), obtained using CosmoNest. 



The prior ranges for the other parameters were chosen 
as in Ref. 0: 0.018 < Q h h 2 < 0.032, 0.04 < tt cdm h 2 < 
0.16, 0.98 < 6 < 1.1, < t < 0.5, and 2.6 < ln(A s x 
10 10 ) < 4.2. Here is a measure of the sound horizon 
at decoupling, and the other symbols have their usual 
meaning. 

When we quote Bayes factors, model Ma is always 
taken to be the Harrison-Zel'dovich case. We normal- 
ize to this case, which means positive numbers indicate 
models preferred against this case. 



A. The spectral index 

For the spectral index ns , we will throughout assume a 
prior range 0.8 < 7i§ < 1.2, as in Refs. [I|l3. The model 
selection results presented for ns must therefore be un- 
derstood in light of this prior. As the allowed regions are 
well contained within this prior, it is trivial to recom- 
pute the Bayes factor if this range is extended; e.g. if it 
is doubled then lni?i is reduced by In 2 ~ 0.7. 

Qualitative understanding of our results can be ob- 
tained from studying the marginalized distributions for 
ns, shown in Fig.^ obtained using CosmoNest. For rea- 
sons explained later, we did not include marginalization 
over the SZ effect in obtaining these likelihoods, which 
shifts them somewhat towards ns = 1 as compared to 
those of the WMAP team Q. A rapid guide to the ex- 
pected result can be obtained by employing a gaussian 
approximation to the marginalized posterior distribution 
for ns. As shown by Trotta [10(, the Bayes factor can 
be computed in this approximation as a function of A, 
being the 'number of sigmas' of the putative detection, 
and the information content / = log 10 (Ans/o') which 
measures the reduction of the allowed parameter volume 
between the prior and posterior (where Ans is the prior 
width and a the standard deviation of the posterior). 
For WMAP3 alone and our choice of prior A ~ 2.3 and 
I ~ 1.4. Employing Eq. (18) of Ref. yields the esti- 
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mate In Bio = 0.4, i.e. the varying ns model is preferred 
but only very mildly. However, this expression assumes 
that the gaussian form holds quite far into its tail, which 
may not be valid, and so we proceed to results from the 
full numerical calculation (which we anyway had to do 
to obtain these marginalized likelihoods as a by-product, 
as explained in Appendix B). 

Using CosmoNest, we ran the Harrison-Zel'dovich 
(HZ) and spectral index cases (ns) to find the difference 
in evidence for WMAP3, with and without the external 
CMB and large-scale structure data sets. We found that 
using WMAP alone, lnB 10 = 0.34, in agreement with 
the estimate above. However, when the extra data sets 
are included, lnJ5io = 1.99, which is substantial but not 
strong evidence for the necessity of ns as an extra pa- 
rameter. The results are given in Table [I] 

The difference of 1.65 in InE between the two datasets 
can be understood simply from the marginalized likeli- 
hoods shown in Fig. ^ Although the curves are similar 
near the peaks, and the maximum likelihood value has 
shifted only by about 0.005, the difference in the mean 
and more importantly the variance have a large effect in 
the tails. The probability of the ns = 1 value is about 4 
times smaller in the WMAP+all case, which would, all 
else being considered equal, translate into a change of 1.4 
in \uE explaining most of the difference. Nevertheless, 
this shift is not particularly significant on the Jeffreys' 
scale. 



B. The inflationary plane 

A non-zero value of ns — 1 is commonly interpreted 
as a strong indication in favour of inflation. However 
slow-roll inflationary models predict not just a non-zero 
ns — 1, but also a non-zero value for the tensor-to-scalar 
ratio r. A proper model comparison motivating inflation 
should therefore examine not the spectral index model 
but the two-parameter extension of HZ into the ns~r 
plane. There is no simple way to make an estimate of 
the outcome in this case. 

At this point we run into the issue of the choice of 
prior for r. While it is uncontroversial to choose a uni- 
form prior for ns, whose value is more or less known, r 
is instead a parameter whose order-of-magnitude is cur- 
rently unknown (sometimes called a 'scale' parameter). 
We will consider two possibilities. The first, Case 1, is 
that r has a uniform prior in the range [0, 1], which is the 
assumption used by Spergel et al. [lj for parameter esti- 
mation. The second, Case 2, considers the Jeffreys' prior 
which states that for scale parameters the prior should 
be uniform in ln(r) rather than r. For the problem to 
be well-defined this needs to be cut off at both ends. 
We use the same upper limit, and as a lower limit take 
the smallest conceivable inflation scale of the electroweak 
scale, which would yield r ~ p 1 ^ 2 /nipj ~ 10~ 34 (where 
p is the energy density). The prior range is therefore 
—80 < lnr < 0. Results are shown in Tabled 



Datasets 


Model 


In Bio 


WMAP only 




HZ 


0.0 






ns 


0.34 ±0.26 


WMAP + all 




HZ 


0.0 






ns 


1.99 ±0.26 




ns + 


r (uniform prior) 


-1.45 ±0.45 




ns 


+ r (log prior) 


1.90 ±0.24 



TABLE I: Evidence differences for the different models and 
different data sets, as discussed in the text. 



For Case 1, CosmoNest calculations indicate Initio 
of —1.45. The large amount of unused prior parameter 
space in the ng-r plane means that this model is some- 
what disfavoured as compared to HZ. 

For Case 2, a calculation is not in fact really neces- 
sary, since the vast majority of the prior space lies in the 
region where r is observationally negligible, and hence 
generates the same likelihood as a model where the spec- 
tral index alone varies. We confirm this explicitly by 
computing evidence over a limited range —5 < lnr < 
and extrapolating the result down to lnr = —80. 

We conclude therefore that the evidence in the infla- 
tionary plane does carry significant prior dependence, 
bracketted by the values we have found under Case 1 
and Case 2. Given the present shape of the likelihood, 
the evidence for the inflation model will not be as large 
as for the spectral index model under any prior choice, 
and may be significantly less. For a uniform prior on 
r, the inflation model is actually rated below Harrison- 
Zel'dovich. 



C. Systematic effects 

The evidence computation we have described takes 
into account only statistical uncertainties. However one 
should also consider the possible effect of systematic un- 
certainties, and there are some indications that these are 
present at a level which would have some impact on our 
conclusions, despite the very careful job that the WMAP 
team have done. We highlight some of these issues here. 

There is some effect from the precise choice of dataset 
used. All the dataset combinations quoted in Ref. 
give very similar constraints on ns, though none corre- 
sponds precisely to the dataset compilation we are em- 
ploying. Curiously though, the dataset WMAP±all on 
the LAMBDA archive, |22j which adds two supernovae 
datasets to our compilation, gives an ns value about one- 
sigma lower than any other dataset quoted, which would 
be expected to lead to a stronger result for the Bayes 
factor. However it is puzzling that this data compilation 
gives a lower ns (and optical depth r) than do any of the 
separate datasets from which it is compiled. 

There is some uncertainty in how to treat the Sunyaev- 
Zel'dovich effect and gravitational lensing. The WMAP 
team allow only for the former, while Lewis has argued 
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22] that the two effects are of the same order, and nearly 
cancel, at least as regards their effect on n§, and that it 
is better to ignore both than to include only one. Ac- 
cordingly we have not included the SZ correction, which 
increases ns as compared to the WMAP3 analysis. 

Another subtlety concerns the modelling of the beam. 
As discussed in Ref. [2l] there are different options for 
doing this, which appear to have a slight effect on the con- 
straint on ng- We have followed the procedure described 
in that paper, rather than that of the main WMAP3 pa- 
pers 0. 

Yet more uncertainty surrounds the modelling of the 
recombination process. According to Ref. HU, inclusion 
of additional two-photon decays leads to significant dif- 
ferences as compared to the standard RECFAST treat- 
ment used in the WMAP papers. If confirmed, this is 
perhaps not too important for WMAP, but would cer- 
tainly matter at Planck sensitivity (Antony Lewis, pri- 
vate communication). 

Also, the reionization optical depth r and ns are cor- 
related. The constraint on r comes mainly from the esti- 
mate of the power in the low I multipoles of CMB polar- 
ization. Substantial foregrounds are present in polariza- 
tion data, so that their removal using just the frequency 
information gathered by WMAP can be tricky. Fore- 
ground subtraction uncertainties could therefore affect r 
and hence ns- 

Finally, we note that the inclusion of Lyman alpha 
power spectrum data (not used in the WMAP3 papers) 
seems to have a marginally significant effect. Accord- 
ing to the analysis of Ref. 24], inclusion of this data 
shifts ns upwards by around one-sigma while leaving the 
uncertainty unchanged. Similar results are obtained in 
Ref. though the trend is less clear as they round 
their quoted results at the second decimal place. 

While individually none of the above would have a very 
major effect on model selection conclusions, that there 
are so many clearly urges caution in interpretting a result 
whose statistical significance remains rather marginal. 



IV. CONCLUSIONS 

We have carried out a Bayesian model selection analy- 
sis of WMAP3 data, as advocated by the WMAP team. 
We have found that WMAP3 data do indeed give support 
for a varying spectral index when combined with other 
data, with the Bayes factor compared to the Harrison- 
Zel'dovich spectrum being approximately Initio = 2. 
According to the Jeffreys' scale, this should be regarded 
as significant, but neither strong nor decisive. It cor- 
responds to probabilistic odds of about 8 to 1 against 
the Harrison-ZePdovich model (i.e. the chance that ns is 
equal to one is about that of tossing a coin three times 
and them all being heads). WMAP3 alone does not pro- 
vide any discrimination between the models. 

In computing our numbers, we have assumed through- 
out that the prior model probabilities are equal, so that 



models are regarded as equally likely before the data 
came along. Anyone who prefers to make an alterna- 
tive assumption is welcome to do so, and can readily 
follow the consequences using the evidence numbers we 
have supplied. For instance, a perfectly plausible stand- 
point might be that since inflation is a physical model, 
its predictions should be taken more seriously than pure 
HZ which is motivated only by symmetry considerations. 
Hence its prior model probability should be greater, per- 
haps tipping the post-data odds decisively against HZ. 
Readers are quite welcome to take that viewpoint, but 
should bear in mind that their conclusion then derives 
from a mixture of the data and their prior prejudice. 
From the data alone, the situation remains to be deci- 
sively resolved. 

In a companion paper j2(j, we forecast the abilities of 
the Planck satellite to resolve this situation, in light of 
the WMAP3 results. 
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APPENDIX A: THE NESTED SAMPLING 
ALGORITHM 

Our implementation of the Nested Sampling algorithm 
is described in Ref. ||. To summarize, the algorithm (as 
first developed in Ref. ) recasts the problem of calcu- 
lating the evidence as a one-dimensional integral in terms 
of the remaining prior mass X, where dX — P(9\M)d6. 
So the integral is transformed 

E = J L(6)p(0)d0 -> E = J L(X)dX (Al) 

where L(X) is the likelihood P(D\9,M). The algorithm 
samples the prior a large number of times, assigning a 
'prior mass' probabilistically to each sample. The sam- 
ples are ordered by likelihood, and the integration follows 
as the sum of the sequence, 

m j 

£ = E^< £i = y(*j-i-*i+i)- (A2) 
j'=i 

The scheme is illustrated in Figure [21 

In order to compute the integral accurately the prior 
mass is logarithmically sampled. We start by randomly 
placing N points uniformly in the prior space, where in a 
typical cosmological application N ~ 300. We then itera- 
tively discard the lowest likelihood point Lj , replacing it 
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FIG. 2: Schematic of the Nested Sampling algorithm. 



with a new point uniformly sampled from the remaining 
prior mass (i.e. with likelihood > Lj). Each time a point 
is discarded the prior mass remaining Xj shrinks by a 
factor that is known probabilistically and the evidence 
is incremented accordingly. In this way the algorithm 
works its way towards the higher likelihood regions. 

As the remaining prior mass shrinks by orders of mag- 
nitude, the challenging part is to find an efficient way to 
draw new points from the remaining prior volume. We 
do this by using the N —1 remaining points at each stage 
to define an ellipsoid that encompasses the extremes of 
the points and is aligned with their principal axes. The 
ellipsoid is expanded by a constant enlargement factor, 
in order to allow for the iso-likclihood contours not being 
exactly elliptical, as well as to take in the edges. New 
points are then selected uniformly within the expanded 
ellipse until one has a likelihood exceeding the old mini- 
mum. 

The process is terminated when the integral has been 
computed to desired accuracy (see Ref. H). In the end 
the evidence contributed by the N — 1 points remaining 
is added to the accumulated evidence. 

The method is general, and the effects of topology and 
dimensionality arc implicitly built into it. 



APPENDIX B: COSMONEST 

The Cosmological Monte Carlo code (CosmoMC) de- 
veloped by Lewis and Bridle j2]j was created to per- 
form an exploration of the cosmological parameter space, 
through the Monte Carlo Markov Cha in p rocess (for an 
overview of MCMC methods see Ref. [3). While it is 
most commonly used with the Metropolis-Hastings algo- 
rithm, other sampling algorithms (such as Gibbs Sam- 
pling and Slice Sampling) can easily be implemented. 



The Nested Sampling algorithm can be considered as just 
another Monte Carlo sampling algorithm. The important 
difference is that the generation of a chain, which in this 
case is not a Markov chain, is ancillary to its primary 
purpose of calculating the evidence accurately. 

The Cosmological Nested Sampling code (CosmoNest) 
we have developed is an additional module that works as 
part of CosmoMC. 



1. Evidence evaluation 

CosmoMC has a 'memory' of only one point: the algo- 
rithm needs only to know where it is in order to decide 
where to go next. CosmoNest needs to know about the 
point it is discarding, but must also hold in its memory 
all the other N — 1 live points, as well as knowing how 
far through the prior mass X it has progressed and what 
value of the evidence (E) it has accumulated. The output 
of a CosmoNest run consists of the set of discarded min- 
imum likelihood points, along with their X value, their 
likelihood, and total accumulated evidence to that point. 

CosmoMC runs multiple chains for two purposes: in- 
creasing the speed of generating samples, and as a way of 
estimating the extent to which the chains have explored 
the parameter space (the Gelman-Rubin statistic) . Here 
we run multiple iterations of CosmoNest to obtain an 
estimate of the uncertainty in the computed evidence. 



2. Posterior samples 

The sequence of discarded points from the Nested Sam- 
pling process is similar to the Markov chain produced by 
an MCMC process with one important difference: the 
MCMC points are sampled from the posterior whereas 
the Nested Sampling points are sampled from the prior 
with a known distribution in X. With the appropriate 
weightings, the 'chain' of discarded points (distributed 
uniformly in \nX) plus the remaining live points (dis- 
tributed uniformly in X within the remaining volume) 
can be used to construct the posterior probab ility distri- 
bution of the parameters, as outlined in Ref. [l4|. 

To summarize, from Bayes' theorem 



p(8\D) = 



L(D\9)ir(d) 
E{D) 



(Bl) 



where p(9\D) is the posterior probability of a parameter 
point 9 given data D, L is the likelihood and ir the prior. 
So for an element i in the chain of discarded points, the 
posterior weighting is 



Pi 



E 



(B2) 



where Wi = ^{X^-y — Xi + i) is the prior mass associ- 
ated with that particular point. The N — 1 points finally 
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FIG. 3: The posterior weights pi assigned to each point in 
one of our HZ runs. The z-axis is the element number in the 
chain, and the vertical dashed line indicates where the live 
points start to be used. This transition could be shifted to 
the right by running the code for longer so as to generate a 
longer chain of discarded points. 

remaining also need to be included to avoid undersam- 
pling the centre of the distribution. They are taken as 
uniformly sampling the remainder of the prior space. 

Figure |3 shows as an example of the posterior weights 
assigned in a particular run. The early points have neg- 
ligible weight as their likelihood is low, and the late ones 
because the prior mass Wi per point becomes small. We 
see that in this case the live points have to be included 
to properly sample the centre of the distribution. The 
fractional contribution from live points can be reduced 
by running the code for longer. The structure of these 
weights should be contrasted with Metropolis-Hastings 
where all samples have integer weights (values greater 
than one accruing when new samples are rejected and 
instead the original sample duplicated). 

Using this method we can reconstruct the posterior 
samples and compare to similar results from standard 
Metropolis-Hastings MCMC. We illustrate this in Fig.[|| 
Posteriors obtained from the two methods are in good 
agreement. 

3. The information 

The information H is defined as (minus) the logarithm 
of the amount the posterior is compressed inside the prior 
[l4|. It is easy to compute from the posterior samples 
once the evidence has been calculated: 

The uncertainty on a single estimate of the evidence is 
dominated by the Poisson uncertainty in the number of 
steps (replacements) to reach the bulk of the posterior. 
This is given by a 2 (E) w H ln[(N + 1)/N}. For the priors 
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FIG. 4: Posterior samples from Nested Sampling (solid) com- 
pared to MCMC (dashed), for a ACDM HZ model using 
WMAP3 data only. 



we have been considering H ~ 10, and given our choice 
of N, this uncertainty turns out to be 0.15 to 0.2. 



4. Public code release 

The CosmoNest code is now freely available for public 
use, and can be downloaded from www.cosmonest.org. 
Its use requires a working installation of the CosmoMC 
package of Lewis and Bridle [27| ■ 
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